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Abstract Using the Wang-Landau flat histogram Monte Carlo (FHMC) simulation technique, we were able 
to study two types of triangulated spherical surface models in which the two-dimensional extrinsic curvature 
energy is assumed in the Hamiltonian. The Gaussian bond potential is also included in the Hamiltonian of 
the first model, but it is replaced by a hard-wall potential in the second model. The results presented in this 
paper are in good agreement with the results previously reported by our group. The transition of surface 
fluctuations and collapsing transition were studied using the canonical Metropolis Monte Carlo simulation 
technique and were found to be of the flrst-order. The results obtained in this paper also show that the 
FHMC technique can be successfully applied to triangulated surface models. It is non-trivial whether the 
technique is applicable or not to surface models because the simulations are performed on relatively large 
surfaces. 

Keywords Triangulated surfaces • Collapsing transition • Surface fluctuations • Flat histogram Monte 
Carlo 



1 Introduction 

Surface models for membranes and strings constructed by Helfrich and Polyakov are described using the 
notion of two-dimensional differential geometry [TJ[2] . The surface shape in is considered to be governed 
by a curvature Hamiltonian, which is given by an integral over the squared mean curvature or the extrinsic 
curvature. Consequently, the surface strength is characterized by bending rigidity b [3l|4l[5l[6l|7] . Thus, we 
understand that the surface collapses and wrinkles in the limit of & — > while it swells and becomes smooth in 
the limit of 6 — ?> cx3 . Theoretical studies utilizing the renormalization group technique predict that the 

crumpling transition is continuous |lllll2l[T31ll4[|15j , while density- matrix renormalization group studies on 
the folding of triangular lattice [16] and recent numerical simulations of the model on triangulated surfaces 
pnfTBlfTOll^ indicate that the crumpling transition is of the first-order and accompanies the collapsing 
transition. 

By including certain inhomogeneous components such as cytoskeletal structure or holes in the above- 
mentioned homogeneous models, we obtain a variety of surface models for numerical studies |21|[^[^ . 
Lateral diffusion of lipids can also be implemented in the models by the so-called dynamical triangulation 
technique, which introduces non- uniform coordination numbers q to the triangulated surfaces [2^. We should 
note that such non-uniform q naturally appears in diagrammatic expansions of the matrix integral in the 
matrix model of 2D quantum gravity, where the surface is embedded in the D = dimensional space, or 
equivalently it is not embedded in any external spaces j25j . 
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In those inhomogeneous surface models in R'^, the transitions separating two neighboring phases are 
discontinuous pTll^f^H] . Thus, the first-order nature of transitions seems to be a common feature of the 
shape transformation transitions in the triangulated surface models. Therefore, careful numerical studies 
are still needed to understand the phase structure of surface models, because first-order transitions are not 
always easy to analyze numerically. 

In this paper, we study conventional homogeneous surface models on triangulated fixed-connectivity 
surfaces of sphere topology using the flat histogram Monte Carlo (FHMC) simulation technique of Wang 
and Landau [5^. The surfaces are allowed to self-intersect. They are called phantom, or self- intersecting, 
surfaces. 

The transitions of homogeneous surface models are already reported to be of the first-order. They were 
obtained using the canonical Metropolis Monte Carlo (MMC) simulation technique, as mentioned above 
[17, 18, 19,20j. MMC is a simple and reliable technique for studying phase transitions in all statistical me- 
chanical systems. However, MMC is not always an efficient technique for analyzing first-order transitions. 
The curvature energy of surface models jumps at the first-order transition point. Consequently, the surface 
configuration can be trapped in one of the minimum energy states if the lattice size increases. The configu- 
rations are trapped because MMC is directly defined based on a canonical ensemble, where the configuration 
of dynamical variables are generated by the Boltzmann weight, which confines the configurations to a narrow 
energy range. 

In contrast to MMC, FHMC is defined by a random walk in energy space. A large number of studies 
have been conducted to verify that FHMC is an efficient and reliable technique for phase transitions in 
spin models [26ll27l[28ll29ll30) . FHMC technique can also be applied to models of polymer chains, in which 
the phase space is non-compact [3T, and also to models with more than two energy terms [321133) . The 
non-compactness of phase space and the multiple energy terms in those models share the common property 
with the surface models. In fact, the phase space of surface models is R^, which is non-compact and is in 
sharp contrast to the compact phase space of spin models. Moreover, the Hamiltonian of surface models is 
always composed of two energy terms: the linear combination of a bond potential Si and a bending energy 
82- However, it is nontrivial whether such a sophisticated technique is useful for studying phase transitions 
in surface models. In fact, a large scale simulation is necessary to study phase transitions of surface models, 
while FHMC technique is considered to be problematic on large systems of a spin model [S^ . 

The purpose of this study is two-fold. The first part aims to determine whether FMHC technique can 
be successfully applied to first-order transitions of surface models. The second purpose is to confirm that 
the conventional, homogeneous models undergo the first-order transition, which was first assessed by MMC 
simulation as mentioned above. 

The density of energy f2 depends on two independent energies, and S2, such that i7 — i7(5i,52). 
Therefore, we must replace ^{81,82) by the single energy density f^{S2), because it seems difficult to obtain 
n(Si,S2) due to the lack of computational speed currently available. Thus, we should check whether this 
replacement is well defined or not. 



2 Models 

The models are defined by the partition function 

Y[dX,exp[^S{X)], (1) 

i=i 

where N is the total number of vertices and S{X) is the Hamiltonian. The symbol YiiLi ^-^i denotes 37V- 
dimensional integration in R'^ . j' indicates that the integrations are performed such that the surface center 
is fixed to remove the translational zero mode. The self-avoiding property of the surface is not assumed. 

Spherical surfaces in R^ are triangulated, and S{X) is defined on them. The triangulated surfaces are 
constructed from the icosahedron by splitting the edges into i pieces of uniform length and dividing the 
faces into triangles accordingly. The surfaces are identical to those used in [TH] . The total number of vertices 
is thus given by 7V= 10£^+2. The total number of bonds Nb and the total number of triangles Nt are given 
by Nb — 3Q£^ and NT — ^Oi"^, respectively. The coordination number q is q = 6 throughout the lattice except 
at 12 vertices, which are the vertices of the icosahedron and of g = 5. 
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The first model denoted by model 1 is defined by the Hamiltonian S{X), which is the linear combination 
of the Gaussian bond potential 5*1 and the bending energy ^2 with the bending rigidity b: 

S{X) = Si+bS2. Si=Y,{X^-X,)\ 

iv) 

^2 -5](l-n, -n,), (modell). (2) 

"^(^ij) in Si denotes the sum over bond (ij) connecting the vertices i and j. j in S2 is the sum over 
triangles i and j, which share a common bond. The symbol {Xi — Xj)'^ in 5*1 is the bond length squares 
between the vertices i and j. The symbol in S2 denotes a unit normal vector of the triangle i. The inner 
product of the normal vectors • rij can also be represented by cos 9ij , where 9ij is the edge angle between 
two triangles. The unit of b is fcT, where k is the Boltzmann constant and T is the temperature. 

The second model denoted by model 2 is defined by the linear combination of a hard wall potential V^g 
and the bending energy S2 with the bending rigidity b such that 

SiX)^Vr„+bS2. 

(ij) 

^2 (model 2), (3) 

(ij) 

where K-o denotes that the potential depends on the parameter tq. The symbol — in Vrg is the 

potential between the vertices i and j and is defined by 

V{\X,-XA) = l^ (|X-X,|<ro), 
^' ^'^ 1 00 (otherwise). ^ ^ 

The value of rp in the right hand side of Eq. ^ is fixed at ro = vTT. Consequently, we have (^{Xi — 
Xj)'^)/N ~ 3/2, which is automatically satisfied in model 1 as described below, where the Gaussian bond 
potential 6*1 = X(Xi — Xj)^ is included in the Hamiltonian in place of the hard- wall potential Vr„. 

The hard wall potential Vrg makes the mean bond length constant just like the Gaussian bond potential 
^i, which makes the mean bond length constant in model 1. If it were not for the constraint jX^— Xj| <ro, the 
size of the surface would grow larger and larger in the MC simulations. Thus the constraint jX^— Xj| <ro is 
necessary to make the bond length well-defined if the Gaussian term Si is not included in the Hamiltonian. 
Consequently, model 2 has an additional parameter, ro, which seems to fix a length scale in the model. 
However, Ref. [31] shows that the results are not dependent on tq. Therefore, we use ro — VTA in the MC 
simulations. 

We should note that (Si/N) =3/2 is satisfied in model 1. This is understood from the scale invariant 
property of Z [TD]. In fact, by rescaling the integration variable in Z such that X — > aX, we obtain 
Z{a) = a^^^~^'^ j' Y[i=i dXi exp [— S'(q;X)], where S{aX) = Si-\-bS2- The scale invariance of Z indicates 
that Z{a) is independent of a and, therefore, is represented by dZ[a) / da\a=i — Q. Thus, we have (Si/N) = 
3(A^— l)/2A^~3/2. In contrast, (Si/N) =3/2 is not satisfied in model 2 because the partition function Z is 
not scale invariant in this case. In fact, the scale transformation X — >■ aX non-trivially changes the potential 
Vrg because tq also changes. Therefore, dZ{a)/da\a=i = is not always satisfied. However, (Si/N) =3/2 is 
almost satisfied in model 2 as we will see below. This finding implies that the violation of the scale invariance 
is not very large in model 2. 

S{X) in Eq. ([2]) naively has an expression S{X) ^ aSi+bS2 with a parameter a of the unit [kT/L'^], 
where [L] is the length unit. By using the above mentioned scale invariance of Z, S{X) can also be written as 
S{X) = Si+bS2- In fact, the partition function Z is unchanged by rescaling X X/ y/a, while this rescaling 
changes S = aSi+bS2 to S ^ Si+bS2- Thus, we can use this simple expression as the Hamiltonian, where we 
should always remind ourselves that the coefficient of Si is assumed to be a[kT/L'^] = 1. 
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3 Flat Histogram Monte Carlo technique 

3.1 Histogram and reweighting 

We first describe the histogram technique for model 1 using the terminologies described in j^S]. As men- 
tioned in the Introduction, the phase space of the models is R^, and the dynamical variables X(g R^) are 
continuous. Consequently, the energies Si and 5*2 of the models are continuous in contrast to the energy 
of the Ising spin model. For this reason, the number of states f]{Si,S2) at energies and ^2 should be 
understood with the symbol dSidS2, such that dSidS2f^{Si, 6*2). Thus, the partition function Zlbo) in Eq. 
dl]) can also be described by i7(5i, 52), such that 



Zibo) ^ JJ dSidS2f2iSi, S2) exp + 60^2)] • (5) 

Let Pbg be Pbg oc f2{Si, S2) exp [—(51 + &0'S'2)]- We then obtain Pf, cx Pb„ exp [— (& — bo)S2]- Using this, we 
have the expectation value of a physical quantity Q at b such that 

{Qm = JJ dSidS2QPb{Si, S2)/ J j dSidS2Pb{Si, S2). (6) 

This expression is a basic formula for reweighting called a single histogram technique, which allows us to 
obtain {Q{b)) from the Monte Carlo data Pbg at bg. 

The reweighting technique described above has two tasks. The first is to obtain f2{Si, S2), and the second 
is to obtain the canonical expectation value by reweighting. The problem lies is performing the first task or 
obtaining f2{Si, S2) efficiently. It is well known that the MMC technique is less efficient for evaluating in a 
large system because of the exponential dumping of f2 in the energy space. The multicanonical Monte Carlo 
simulation (MCMC) technique is considered a dynamic version of the multihistogram reweighting technique 
|26| . However, it is not apparent whether this technique is applicable to the simulations of large surfaces. 
In fact, our preliminary study indicates that a flat histogram is barely obtained by MCMC simulations on 
the surfaces of A^>5762 for model 1. By contrast, a fiat histogram can be obtained even on relatively large 
surfaces by the fiat-histogram MC technique, which was recently proposed by Wang and Landau. Detailed 
information on the FHMC technique for surface simulation is described in the following subsection. 



3.2 Flat histogram Monte Carlo 

The Wang-Landau FHMC technique for surface models, which is simply denoted as FHMC, consists of the 
following three steps [26] : 

1) Recursive construction of ]7(52). 

2) A production run with J7(52), collecting measurements. 

3) Reweighting to extract the canonical expectation of physical quantities. 

The first and the second steps are long Monte Carlo runs, but the third step is very short by comparison. 
The second run (=Step 2) includes thermalization Monte Carlo sweeps (MCS). In the following, we describe 
each step in detail. 

In FHMC simulations for model 1 and model 2, we use a reduced density of states ^{82), which is 
obtained by the replacement 

0(81,82)^ ^{82). (7) 

This replacement is necessary because ^2{8l,82) is a double histogram and, hence, seems very hard to obtain 
accurately. This replacement is also well defined if we recall that 81 is almost constant in model 1 due to 
the scale invariance of the partition function. The fact that 5*1 is almost constant implies that the role of 
81 is to make the mean bond length constant in model 1. Therefore, the potential 5*1 can be replaced by a 
Lennard- Jones type potential |17j or by a hard-wall potential, which is the potential in model 2. Thus, the 
phase structure of the surface models is primarily dependent on the curvature energy 5*2 . For this reason, 
the replacement in Eq. ([7]) is considered as reasonable. The histogram 0(82) is defined by 82 in the region 
5'™'" < 5*2 < 8^^^, where 5™'" and S'™^'^ are chosen such that the phase transition region is included. We 
should note that ^2(82) implicitly depends on 81 in model 1 and on Vrg in model 2, although the replacement 
of Eq. d?]) is assumed. 
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In Step 1, ^2(82) is recursively obtained by updating the variables X as follows: the new position of 
the vertex i is given by X^~Xi + AX, where AX is chosen randomly in a small sphere. The new position 
X- is accepted with the probability Min[l, Q{S2)/ f2{S2)], where S2 and 5*2 are given by S2 = 6*2 (old) and 
^2 = 5*2 (new). The radius of the small sphere is chosen so that the acceptance rate rx for X' is approximately 
rx =50%. However, rx varies during the simulations in Steps 1 and 2 in contrast to the MMC case, where 
rx remains almost constant. The constraint 3/2— Zi<S'i(new)/A^<3/2 is also imposed on X^ only in model 
1, while the constraint V in Eq. Q is imposed on X^ in model 2. At the beginning of the simulation of Step 
1, MMC simulations are performed as a thermalization MCS to make the variables X satisfy the conditions 
3/2- A < Si (new) /N<3/2 and S'^"' < 5*2 < S'i"^''. The parameter A is fixed to Zi = 0.05 in model 1. This 
constraint for 6*1 (new) /N is imposed on the new variable X^ so that the replacement of Eq. (O becomes 
well defined, as mentioned above. 

The histogram 17(5*2) is updated such that 

n{S2) f n{S2) (8) 

with the initial values 0(82) = 1 and / = exp(l). The update of f2{S2) in Eq. ^ is performed independently, 
whether X' is accepted or not, in every update of X. The formula log[i7(S'2)] — ?> log(/) + log[J7(S'2)] is used 
in the simulations. The multiplicative factor / in Eq. ([S]) is also redefined such that 

/ ^ v7 (9) 

only when a given condition is satisfied. To perform this redefinition of /, we check whether the energy 
histogram H{S2) is sufficiently flat or not at every 10^ MCS by the condition H{S2) < eH, where H is the 
mean value oi H{S2) and e = 0.9. Once this flatness condition for -ff (^2) is satisfied, the histogram H{S2) is 
reset to zero for all 5*2. This procedure is repeated while / > 1 + 10~* in Step 1. The total number of MCSi 
for the recursion is not an input parameter, but an output datum, and its value depends mainly on the 
conditions of / and H{S2)- The parameters 5™™, S^'^^, and As should be provided as input data for the 
recursion simulations so that the histograms f2{S2) and H{S2) are defined at Sf"\ Sf"^ + As, Sf"^+ 2As, 
• • •, S^'^^ — As. Let Nh be the total number of these energy points. Then, the energy step As is given by 
As = {Sf''''-S^''')/NH. 

Step 2 is the production run, which is performed using the histogram f2[S2) and is simply described as 
follows: the new position X[ of the vertex i is accepted with the probability Min[l, exp(— ZiS'i)f2(S'2)/i?(<S'2)], 
where zi^i = 5i(new)-S'i(old), and 12(52) = J7[S'2(old)], /?(5^) = /?[52(new)]. The factor exp(-zi5i) is set 
to 1 in model 2. Because 17(52) is defined in the region 5"'"<52 <S™^^, the new position X[ is limited so 
that 5™™ < 52 < 5™'*'^ just like in Step 1. In Step 2, no constraint is imposed on the bond length in model 

1, while the potential V in Eq. ([3]) imposes its constraint on the bond length in model 2. In both models, 
5i is expected to be almost constant, such that Si/N = 3/2. The thermalization MCS is fixed to 2 x 10^ ~ 
5 X 10^ in both models. A sufficiently large number of MCS2 for production runs is performed after the 
thermalization MCS, where MCS2 is an input parameter in contrast to MCSi in Step 1. The measurements 
are performed every 500 MCS in both models. The acceptance rate rx varies during the simulations in Step 

2, just like in Step 1 mentioned above. 

Table[T]shows the parameters, including 5™'" and S^'^^. Nh is the total number of energy points for the 
histograms of 17(52) and H{S2). MCSi and MCS2 are the total number of MCSs performed in Step 1 and 
Step 2, respectively. 

Step 3 is performed by using a technique that is analogous to the multi-histogram reweighting technique; 
however, the recursion is not necessary because of the single Monte Carlo simulation performed in Step 2. 
The canonical expectation value of a physical quantity Q(b) is given by 

E5. Q hiS2, QMS2) exp(-b52) 
^^'^ Es.EQM^2,Q)r2(52)exp(-652) ' ^^"^ 

where h{S2,Q) is the histogram of the data obtained in Step 2. This formula Eq. ^TU\i corresponds to Eq. 

where 60 is actually set to zero. We should note that the sum over Si in the denominator/numerator is 
dropped from Eq. ((TU)) in model 1. This is because of the replacement of Eq. ([7]). A normalization constant 
can be included in the exponential factor in both the denominator and the numerator in Eq. (jlOp . 

We use a random number called Mersenne Twister [35) . A sequence of uniform random numbers is used 
for a 3-dimensional move of the vertices X and for FHMC accept/reject decisions in the update of X. 



6 



Table 1 

MCSi. 



The parameters used in the FHMC of model 1 and model 2. These are all input parameters, excluding 



model 


iV 






N„ 


MCSi(xlO*) 


MCS2(xlO*) 


1 


2562 


0.35 


0.47 


1000 


1.8 


9 


1 


5762 


0.35 


0.5 


2500 


11.2 


17.8 


1 


7292 


0.35 


0.5 


3300 


15.3 


27.3 


1 


10242 


0.35 


0.5 


4500 


18.6 


49 


1 


15212 


0.35 


0.5 


6700 


44.5 


50 


2 


2562 


0.34 


0.49 


1000 


6.6 


9 


2 


4842 


0.35 


0.49 


1500 


16.8 


21.6 


2 


8412 


0.35 


0.51 


2500 


43.2 


37.2 


2 


16812 


0.35 


0.51 


5550 


45.3 


108 



4 Results 

The mean square size is defined by 



X 



N ^ ' 



(11) 



where denotes the sum over vertices i, and X is the center of the surface. X"^ is identical to the radius 
squares if the surface becomes smooth and spherical, while it becomes X^ if the surface collapses. 
Therefore, X^ can reflect the surface size. 



X 



100 



model 1 



N=15212 

; N= 10242 




N=5762 



N=2562 




0.77 , 0.78 

b 




0.77 , 0.78 

b 



N=2562 



0.685 



0.69 , 0.695 

b 



■ + ' 


. (d) 


model 2 






[ N=8412 








pN=4842 - 








N=2562 ■ 

V : 


N=16812 ^ 






0.685 


0.69 


, 0.695 

b 



Fig. 1 (a) The mean square size vs. h of model 1, (b) the bending energy S2/NB vs. b of model 1, (c) vs. 6 
of model 2, and (d) S2/N B vs. b of model 2. The solid curves are the results of FHMC, while the symbols (Q, 
• • ■) are those of MMC in IT9l[20l. 



Figures [Ija) and [T]Jb) show X^ vs. b and the bending energy S2/NB vs. b of model 1, respectively. 
Nb{=3N — 6) is the total number of bonds. The solid curves are the results of FHMC, while the symbols 
such as the circle, triangle, etc., are those of MMC, which are identical to the results reported in [TO]. The 
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error bars for the symbols in this paper are standard errors, which are obtained by binning analysis. From 
these symbols, we can directly see whether the results of FHMC coincide with those of MMC. For this 
reason, MMC data from [T^ are presented in this paper. Figures [TJc) and[TJd) show the results of model 2. 
The symbols are the results of MMC reported in [5D]. In the case of model 1, we see that the FHMC results 
are in good agreement with MMC in both and S2/NB at iV < 10242. On the TV = 15212 surface, the 
results are only slightly different from each other. The FHMC results of model 2 are also in good agreement 
with those of MMC at iV < 8412. A deviation can be seen in the data obtained on the largest surface of 
TV = 16812. We have no definite explanation for these deviations. However, the total number of MMCs was 
not sufficiently large, at least on the largest surfaces in model 1 [19] and in model 2 [20]. Note also that 
the FHMC results seem to be dependent on the choice of S™™ /Nb and S™^^/Nb if the range is narrow. A 
wide range is better for accuracy in the obtained results. However, the simulation time becomes longer and 
longer with increasing range on such large surfaces. 




Fig. 2 (a) The variance C^a vs. b of model 1, (b) the specific heat Csj vs. b of model 1, (c) C_y2 vs. b of model 2, 
and (d) Csj vs. b of model 2. The solid curves are the results obtained by FHMC, while the symbols are those of 
MMC. 



The variance Cx^ of X'^ is defined by 

Cx^^^{{X^-{X^)f). (12) 

This variance Cx^ refiects the fiuctuations of X^ around the mean values. The specific heat of 5*2 is 
defined by 

Cs.^p{S2~{S2)f). (13) 

This also refiects the fluctuations of the bending energy. 

Figure dja) shows Cx^ vs. b of model 1, while Fig. [Hb) shows vs. b of model 1. The results of 
model 2 are also shown in Figs, djc) andjU^d). The solid curves in the figures are the results obtained by 
FHMC, while the symbols are those obtained by MMC, as presented in Figs. [Ha)-md). In Fig. [2ljb), the 
solid curves of C52 almost agree with the symbols. However, there is a small difference between the two 
results at = 15212. The specific heat obtained by FHMC is not always in good agreement with those 
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obtained by MMC, as can be seen in model 2 in Fig.[2l^d). We feel that the reason for this deviation is mainly 
due to a lack of statistics in MMC simulations in [TOll^ . The MMC data of C^^ in [Ml^ are not shown 
in Figs. [5Ja) and [2Jc) because the deviation between the MMC data and the solid lines is relatively large 
compared to those in Csj. The variance Cx^ and the specific heat are not so easy to obtain accurately. 
However, we observe that the peaks of Cx^ and grow with increasing TV. Therefore, our expectation is 
that the phase transitions are reconfirmed by FHMC simulations. 

Cmax 
X2 



0.05 

0.01 
0.005 



1000 iv 500010000 1000 KT 500010000 

Cmax ^ y^max 

X2 

0.1 
0.01 



500010000 N 500010000 

Fig. 3 Log-log plots of the peak values C^f vs. A'^ and Cg^'"' vs. iV of model 1 and model 2. The error bars for the 
data are standar d e rrors. The data and the errors were obtained by FHMC. The straight lines are drawn by fitting 
the data to Eq. ()14[). The largest three data points are used in the fitting in the cases of (b) and (d). 

The peak values C^f^ of model 1 are shown as symbols in Fig. [3ja) against N in a log-log scale. The 
peaks Cf^^ of model 1 are shown in Fig. [3fb) . The error bars drawn on the symbols denote the standard 
errors of C^f^ and Cg^^. The straight lines in Figs.[3][a) andEJb) are drawn by fitting the data to 

C^^^'^N", Cf^^'^N^, (14) 

where a and /i are critical exponents. The fitting of C^2^ in Fig. [Sljb) is performed by using the largest 
three data points. Figures [31Jc) and|21[d) are the results obtained by model 2. The symbols and the errors 
are simulation data. The straight lines are drawn by fitting the data. The fitted line in Fig. EJd) is obtained 
using the largest three data points. Thus we have 

(T = 1.24 ±0.07, = 0.98 ±0.12, (modell), 

(7 = 1.42 ±0.02, ^=1.26 ±0.05, (model 2). (15) 

The exponent ji of model 1 is almost identical to = 1. Therefore, the transition of surface fluctuations 
in model 1 is considered to be of the first order. We also see that the remaining exponents are consistent 
with a first-order collapsing transition and a first-order transition of surface fluctuations, although they are 
slightly larger than 1. These results are consistent with the previous conclusions described in We 
should emphasize that the symbols and the errors in Figs. [3ja)-l3l[d) are the ones obtained by the FHMC 
technique in contrast to those shown in the previous figures, such as Fig. [T] and Fig. [21 where the symbols 
correspond to MMC simulation data in P^I20| . 
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Si/N 
1.502 

1.5 

1.498 



(a) 



model 1 



0:N=15212 
A:N=10242 
□ :N=5762 




0:N=16812 
A:N=8412 
□ :N=4842 
v:N=2562 



0.765 0.77 , 0.775 

b 



0.685 0.69, 0.695 
b 



Fig. 4 The Gaussian bond potential Si/N vs. b of (a) model 1 and (b) model 2. The solid curves are obtained by 
FHMC, while the symbols are MMC data. 



Finally in this subsection, we show the Gaussian bond potential Si/N vs. b in Figs. Sl^a) andlD^b). The 
solid curves represent FHMC results, while the symbols are the results of MMC, just as presented in Figs. 
[T]and[21 As described in Section^ Si/N in model 1 is expected to be Si/N ^3/2. We see from Fig. HJ^a) 
that this expectation is satisfied just like the MMC data. The solid curves of Si/N of model 2 in Fig. UJ^b) 
are also consistent to the results of MMC. We should note that the discontinuity seen in Si/N reflects a 
discontinuous change of bond length in model 2. Therefore, the phase transition in model 2 is considered to 
be accompanied by a structural change of its surface, although the discontinuity is very small compared to 
the value of Si/N. The fact that the discontinuity in Si/N is very small compared to Si/N itself indicates 
that the violation of scale invariance is small in model 2, as mentioned in Section 2. 



5 Summary and conclusions 

In this paper, we studied two types of triangulated spherical surface models using the flat histogram Monte 
Carlo (FHMC) simulation technique of Wang and Landau. The surface models have long been studied 
numerically, mainly using Metropolis Monte Carlo (MMC) simulations. It has been recently reported that 
the surface fluctuations in these models undergo a first-order transition. We aimed to confirm that the 
transition is of the first order and to verify that the FHMC technique can be applied to study phase 
transitions of surface models. 

It is nontrivial whether FHMC can be used to study surface models because large-scale simulations are 
necessary to study the phase structure of surfaces. It was pointed out in [S^ that the FHMC technique 
can not always be successfully applied to large-scale simulations of spin models. On the other hand, the 
variables X of the surface models are updated by a random walk in energy space in the FHMC technique. 
Therefore, the FHMC technique is considered to have an advantage in studying first-order phase transitions, 
as compared to the MMC technique. However, the surface simulations by FHMC remained to be studied. 

The FHMC simulations performed in this paper consist of three steps: the first is a random walk in 
energy space to obtain the density of states ^7(52) recursively. The second is to collect measurements using 
f2{S2)- The final step is to obtain canonical expectations using a reweighting technique. The first two steps 
are long MC runs, each of which is a single MC run. We need many parameters to start FHMC simulations. 
Short MMC simulations are required to obtain these parameters before starting FHMC simulations. 

The results obtained by FHMC simulations were compared to those reported in [T?ll20j . which were 
obtained by MMC simulations. FHMC results were found to be in good agreement with MMC results in 
both model 1 and model 2, excluding the data obtained for the largest surfaces in both models. Although the 
variances Cx^ and Cs^ obtained by FHMC simulations are slightly different from those previously obtained 
by MMC simulations on the largest surfaces, the order of the transitions is not infiuenced. The reason for 
these differences seems to be due to the lack of statistics in the previous MMC simulations, at least on 
the largest surfaces. No problem is observed in FHMC simulations concerning the system size at least on 
the surfaces upto N — 15000 ^ 17000. Thus, we conclude that the first-order nature of the transitions 
is reconfirmed by the FHMC technique and, consequently, that the FHMC technique can be successfully 
applied to study the phase structure of surface models. 
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